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ANALYSIS  OF  LIFE  TABLES  WITH  GROUPING  AND  WITHDRAWALS 


by 

Dennis  V.  Lindley 

We  use  the  model  of  Breslow  and  Crowley  (1974).  Associated  with 
each  individual  is  a pair  of  independent  random  quantities,  X , the 
lifetime  and  Y , the  withdrawal  time.  The  raw  observations  for  an 
individual  are  Z = min  (X,Y)  , the  time  at  which  he  leaves  either 
through  death  or  withdrawal,  and  an  indicator  which  says  whether  the 
departure  was  caused  by  death  or  withdrawal.  The  time  scale  is  then 
divided  into  nonoverlapping  intervals  and  the  quantity  Z grouped  so 
that  observation  is  only  made  on  the  interval  within  which  he  left  the 
system.  The  quantities  for  different  individuals  are  independent  and 
identically  distributed.  Consequently  if  N individuals  are  present  at 
the  beginning  of  an  interval,  the  data  consists  of  D , the  number  who 
were  observed  to  die  in  the  interval;  W the  number  observed  to  withdraw 
alive  during  the  interval;  and  S the  number  who  survived  to  enter  the 
next  interval.  N=D+W+S.  In  this  paper  we  shall  only  consider  a 
single  interval,  except  for  a remark  in  the  discussion  at  the  end  of  the 
paper.  Our  task  is  to  estimate  the  chance  of  death  occurring  in  the 
interval.  If  X and  Y have  distribution  functions  F and  H respective- 
ly and  the  interval  is  (0,1]  conventionally,  then  this  chance  is  F(l)  . 

The  novelty  lies  in  the  fact  that  the  probability  of  being  observed  to  die 
is  necessarily  less  than  F(l)  , since  some  deaths  during  the  interval  will 
take  place  after  withdrawal  and  therefore  be  unobserved.  This  happens 
whenever  Y < X 1 . The  unknown  quantity  of  interest  is  F(l)  : the 
data  are  D , W and  S , so  that  we  require  to  calculate  p(F(l)  | D,W,S)  . 
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Another  possible  interpretation  of  the  model  is  to  replace  "death"  by 
"death  from  a specific  cause,"  say  cancer;  and  "withdrawal"  by 
"death  from  other  causes."  In  that  situation  we  have  a competing  risk 
problem  and  both  death-rates,  F(l)  and  H(l)  , will  be  of  interest. 
l.e  shall  therefore  concentrate  on  the  calculation  of  p(F(l) ,H(1)  | D,W,S) 
from  which  the  earlier  quantity  can  be  obtained  by  taking  the  marginal 
distrubution.  Notice  that  this  competing-risk  model  assumes  the  two 
risks,  expressed  through  X and  Y , are  independent.  The  model  has 
applications  to  life-tables  in  biology  and  in  reliability  engineering. 

The  model  is  specified  in  terms  of  F(x)  and  H(x)  for  0 < x < 1 , 
with  F(0)  = H(0)  = 0 , so  that  it  is  natural  to  extend  the  conversation 
from  the  quantities  of  interest,  F(l)  and  H(l)  , to  include  these 
functions.  We  therefore  begin  by  calculating  the  likelihood  of  F(x)  , 
H(x)  , 0 ^_x  < 1 for  the  data  (D,W,S)  . An  individual  will  be  observed 
to  die  in  the  interval  if  he  dies  at  any  time  x , 0 < x <_  1 , and  his 
withdrawal  quantity  Y exceeds  x . Hence 


(1) 


p(observed  death) 


[1  - H(x)]dF(x) 


* F(l)  - 


/ 


H(x)dF(x) 


0 


demonstrating  that  the  true  death  rate  exceeds  the  chance  of  observed 
death  by  a nonnegative  integral.  Similarly 
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p(observed  withdrawal)  = J [1  - F(x)]dH(x) 

0 

(2)  i 

» H(l)  - J F(x)dH(x)  . 
0 


To  survive  he  must  neither  die  nor  withdraw,  hence 


(3) 


p(survival)  = [1  - F(l)][l  - H(l)]  . 


It  is  easy  to  verify  that  these  probabilities  add  to  one,  (This  provides 
an  unusual  proof  of  the  usual  formula  for  integration  by  parts!) 

To  economize  on  notation  we  write  F(l)  = <$>  » H(l)  = 9 and  define 


(4) 


P = J F(x)dH(x)/F(l)H(l) 
0 


A probabilistic  interpretation  of  p is  that  it  is  the  conditional 
probability  of  death,  given  that  both  X and  Y are  less  than  one: 
that  is,  given  that  both  withdrawal  and  death  take  place  in  the  interval. 
It  is  an  artifice  of  the  model,  rather  than  a physically  interpretable 
quantity.  Similarly  (1  - p)  is  the  probability  of  withdrawal  under 
the  same  conditions.  Then 


(5) 


r "" ' 


p (observed  withdrawal)  * 0(1  - P<|>)  , and 
p(survival)  = (1  — 4>)(1  — 9)  . 


The  likelihood  for  the  data  is  consequently 


(6)  *D(1  - p$)W(l  - <fr)S{l  - (1  - p)9}D0W(l  - 0)S  . 


Within  the  restrictions  of  the  model  this  expression  contains  all  the 
information  in  the  data  and  several  important  remarks  can  be  made  about 
it.  First,  although  the  conversation  was  extended  from  F(l)  and  H(l) 
to  include  these  functions  throughout  the  whole  of  the  interval,  in 
fact  the  only  additional  quantity  thereby  introduced  is  p . Consequently 
the  data  only  contains  information  on  three  aspects  of  F and  H , 
namely  <}>  , 9 and  p . The  value  of  p will  be  considered  below. 

Second,  the  trinomial  likelihood  (6)  depends  on  these  three  parameters 
<f>  , 0 and  p , each  of  which,  irrespective  of  the  values  of  the  others, 
can  assume  any  value  in  the  unit  interval.  Consequently  in  the  unit 
cube  of  (<)> , 0 , p ) -values , there  is  a one-dimensional  curve  on  which,  as 
<)>  , 0 and  p vary,  the  probabilities  (5)  remain  unaltered,  and  we  have 
a problem  that  econometricians  describe  as  unidentified:  that  is,  even 
if  the  probabilities  on  the  left-hand  sides  of  (5)  were  known  exactly 
(N  -+  ®)  the  values  of  0 , <p  and  p would  not  be  known  but  only  known 
to  lie  on  the  curve  just  mentioned.  Third,  if  p were  known,  then  the 
likelihood  (6),  now  dependent  only  on  9 and  <p  , factors  into  a product 
of  two  terms  which  depend  on  0 alone  auu  <f>  alone.  Consequently 
estimation  of  each  can  proceed  separately  if  9 and  <j>  are  judged 
independent.  In  particular,  as  some  writers  suppose,  if  the  withdrawal 
times  are  deterministic  then  only  the  cfi-part  is  relevant.  But  with  p 
unknown,  the  situation  is  more  complicated  and  there  is  no  obvious  reason 
for  discarding  the  0-component,  when  intersected  in  <(>  , since  it 


appears  to  contain  information  about  p which  may  be  of  value  in  inter- 
preting the  ^-component . Hence  p , which  arises  from  the  grouping, 
plays  a key  role  in  the  analysis,  and  we  now  discuss  it. 


5 


Since  F and  H are  both  increasing,  0 < p £ 1 . Both  extremes 
are  attainable:  for  if  all  withdrawals  take  place  at  the  beginning  (end) 
of  the  interval,  but  no  deaths  occur  there,  then  p * 0(1)  . If  linear 
splines  are  used  to  approximate  F and  H , so  that  *-fey  are  both  linear 
in  (0,1]  , then  simple  calculation  shows  p * \ . More  generally, 
if  F(x)  **  cH(x)  , then  p = h , for  any  c : in  particular  if  c * 1 
and  they  are  identical.  If  deaths  and  withdrawals  take  place  at  constant 
rates,  X and  y respectively,  during  the  interval,  so  that  F(x)  * 1 - e ^ 

_ 1i  V 

and  H(x)  = 1 - e , then  simple  calculation  shows  that 

(7)  „ . > - <>  + “>«~y  . 

(X  + u)(l  - e~*)(l  - e"P) 

Expanding  both  numerator  and  denominator  in  powers  of  X and  y up  to  and 
including  third  powers  shows  that  p * % for  all  X and  y . Hence 
unless  the  death-  and  withdrawal-rates  are  high,  p - h is  a good  approxi- 
mation. Other  assumptions  are  possible:  for  example,  Chiang  (1961)  supposes 
that  all  withdrawals  take  place  at  the  mid-point  of  the  interval  and  that 
the  deaths  occur  at  a constant  rate.  Then  p = F(J5)/F(1)  = (1  - (1  - 

These  considerations,  allied  to  the  unidentif iability  aspect,  mean 
that  two  types  of  procedure  are  possible.  First,  assume  enough  about  F 
and  H for  p either  to  be  known  or  to  be  known  as  a function  of  <p  and 
6 . (Thus  linearity  gives  p **  h : Chiang' s procedure  is  an  example  of 
the  second  possibility.)  Second,  leave  F and  H sufficiently  general 
for  p to  be  unknown,  even  as  functions  of  p and  0 , and  to  face  up  to 
the  identif iability  problem.  Previous  work  on  the  problem  has  concentrated 
on  the  first  procedure,  making  assumptions  about  F and  G . However 
Breslow  and  Crowley  (1974)  show  that  if  this  is  done  then  the  resulting 
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estimators  of  <j>  and  0 are  typically  inconsistent,  in  the  sense 
that  as  the  size  of  the  data  base  increases,  N •+•  °°  , they  do  not  tend 
to  the  true  death-  and  withdrawal-rates.  For  example,  the  classical 
estimate,  D/(N  - %W)  is  only  consistent  under  rather  special  assumptions 
about  F and  H . In  this  paper  we  therefore  take  the  second  view, 
that  p is  not  known,  even  as  a function  of  $ and  0 . We  show  that 
estimation  is  still  possible  and  that  the  possible  inconsistency  is 
avoided . 

We  therefore  return  to  the  likelihood,  given  by  (6),  in  which  , 0 
and  p can  take  all  values  in  the  unit  cube.  Since  p is  unknown  it 
has  a distribution,  which  we  suppose  independent  of  $ and  0 . Let  its 
density  be  f(p)  . (These  assumptions  will  be  discussed  below  when  we  are 
in  a better  position  to  appreciate  the  roles  they  play.)  Consequently, 
from  (6) 

1 

(8)  p(D,W,S  j 0,<fr)  = *D(1  - <P)1 * * * SeW(l  - 0)S  J (1  - p0)W{l  - (1  - p)0}Df(p)dp  . 

0 

Even  with  simple  forms  for  f(p)  , for  example  Beta,  the  integral  cannot 
be  evaluated  in  terms  of  standard  functions.  There  are  two  possible 
procedures:  for  numerical  values  of  D , S and  W to  use  numerical 
integration  techniques;  for  large  N to  use  analytic  approximations. 

In  the  present  paper  we  explore  the  latter  possibility.  This  will  give 
us  the  "feel"  of  the  situation  more  easily  than  the  computational  method, 
though  that  method  will  be  essential  for  small  samples. 

If  N , and  consequently  D and  W , are  large,  the  integrand 
in  (8) , as  a function  of  p , has  a maximum  where 


1 
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or 

(9) 


1 - (1  - P)9  _ 1 - p» 

D0  W<(> 


P 


d _ w(l  - 6) 
0 


n 


say 


| 


where  d = D/(D  + W)  and  w = W/ (D  + W)  ; d + w = 1 . (The  interpretation 
of  d is  the  proportion  of  those  seen  to  die  out  of  all  who  left  the 
system  in  the  interval.)  This  value  of  p , that  has  been  written 
n = n(4>>9)  , is  the  maximum  likelihood  value  of  p for  given  $ and  0 . 

There  are  two  possibilities:  if  n lies  outside  the  range  of 
integration  in  (8),  namely  the  unit  intervals  then  the  integrand  will  be 
small;  if  0 < n < 1 then  the  integrand  will  be  appreciable  only  in 
the  neighborhood  of  the  maximum.  We  shall  therefore  suppose  the  likelihood 
vanishes  outside  0 < n < 1 . 

We  next  use  the  result  that  if  0 < n < 1 

1 

(10)  J eNL(p)f(p)dp  ~ eNL(Tl)f(n){-L"(n)}"1'* 

0 

asymptotically  as  N 00  . (The  double  prime  denotes  a second  derivative.) 
This  is  easily  established  by  a s teepest-descent  type  of  argument. 

On  applying  this  to  «'8)  with 

NL(p)  = W log  (1  - p<j>)  + D log  {1  - (1  - p ) 0 } 

and  simplifying  the  result,  we  have  for  p(D,W,S  | 0,<J>) 


(1  - $)s(l  - 0)s[l  - (1  - e ) (l  - ^)]N"S+1f(n)0_1^"1 
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A change  of  variables  from  $ and  0 to  C = (1  - 0)(1  - $)  and  n 
gives 

(li)  p (d , w, s I c,n)  = ;S(1  - 0N_Sf(n)  • (l  - c)/e<t>  . 

It  may  help  in  the  appreciation  of  what  is  happening  to  consider 
p(D,W,S  | 5,n)  in  the  plane  of  $ and  0 (see  the  figure).  It  is 
zero  (or  strictly,  very  small)  except  between  the  two  curves  n = 0 and 
n = 1 drawn  in  the  figure.  Between  these  two  curves  the  function  is 
o ily  large  in  the  neighborhood  of  the  curve  £ = S/N  , also  drawn  in 
the  figure.  Consequently  we  need  only  concern  ourselves  with  this 
probability  in  a region  indicated  by  hatching  in  the  figure.  Within  this 
region  and  along  curves  £ = constant,  the  probability  is  dominated  by 
f (n)  where  f is  the  prior  density  for  p . Furthermore  if  the  final 
term  in  (11),  (1  - C ) / © > which  comes  from  the  second  derivative  in  (10), 

is  ignored,  the  probability  (or  more  correctly,  likelihood)  factors  into 
functions  of  £ and  n separately.  The  curve  X,  = S/N  intersects  the 
two  curves  n = 0 and  n = 1 at  A and  B in  the  figure,  whose 
^-coordinates  are  respectively  D/N  and  D / (N  - W)  . Hence  with  probability 
almost  one,  $ lies  between  these  two  values  which  corresponds  to  assuming 
that  all  withdrawals  take  place  at  the  end  (beginning)  of  the  interval. 

A common  estimate  is  to  adopt  the  compromise  D/ (N  - '^W)  . Most  of  these 
results  can  be  deduced  by  equating  the  chances  in  (5)  to  D/N  , W/N  and 
S/N  respectively.  The  additional  feature  that  our  analysis  shows  is  the 
factorization  of  the  likelihood  when  expressed  in  terms  of  z,  and  n , 
and  the  explicit  form  of  the  factors,  in  particular  that  in  n , which  is 
entirely  due  to  the  prior  opinions  about  p , f(*)  . 


If,  over  the  hatched  region,  the  prior  density  for  4>  and  9 is 
reasonably  constant  then  the  expression  on  the  right-hand  side  of  (11) 
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is  approximately  proportional  to  the  probability  of  £ and  n , given 
D , W and  S . If  the  contribution  from  the  second  derivative  is 
ignored,  we  have  approximately 


p(C,n  | D,W,S)  - £S(1  - c)N_sf (n) 


so  that  (;  and  n are  independent,  C,  having  a Beta-distribution 
and  n having  the  density  of  p . Admittedly  rather  crude  approximations 
have  been  made  in  order  to  obtain  this  result,  but  experience  with  steepest- 
descent  type  arguments  in  this  sort  of  situation  suggests  that  they  will 
work  reasonably  well  and  not  demand  excessi/ely  large  data  sets  (large 
values  of  N) . They  are  about  as  good  as  maximum  likelihood  results  are. 

The  posterior  mode  of  (12)  occurs  at  C = s and  n = r , where 
s = S/N  , the  proportion  who  survive  the  interval,  and  r is  the  most 
probable  a priori  value  of  p . Translating  this  back  into  values  of 
<j>  and  9 , the  modes  satisfy 


(1-40(1-9)  = s and 


d w(l  - 9) 


<{)  9 


Solving  for  41  > the  chance  of  death  in  the  interval,  we  easily  obtain 
the  quadratic  equation 


"r  - 4>{  (d  + r)  + s(w  - r)}  + d(l  - s)  = 0 . 


In  terms  of  the  raw  data  D , W and  N , this  is 
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(15)  <J>2 * *rN  - <J>{N  - W + r(D  + W) } + D = 0 . 

It  was  remarked  in  the  third  comment  after  (6)  that  if  p is  known 
the  <})-part  of  the  likelihood  is 


*DU  - P*)Wd  - *)S  • 

It  is  easy  to  verify  that  the  maximum  of  this  likelihood  satisfies  (15) 
with  P = r . In  other  words,  our  posterior  mode  for  <j>  is  the  maximum 
of  the  <J>-part  of  the  likelihood  with  p replaced  by  its  most  likely 
value,  r . The  left-hand  side  of  (14)  equals  d(l  - s)  > 0 at  <j>  = 0 
and  equals  -s(l  - r)  < 0 at  <t>  = 1 , so  that  there  exists  just  one 
root  in  the  unit  interval.  At  the  value  of  <)>  given  by 


(16) 


D 

N - (1  - r)W 


2 3 

the  left-hand  side  of  (15)  is  easily  calculated  to  be  r(l  - r)D  W/N  , 
which  is  small.  Hence  in  many  cases  is  a reasonable  estimate  of 

the  death-rate.  It  can  be  justified  on  intuitive  grounds  by  supposing 
that  all  withdrawals  take  place  at  r , so  that  each  withdrawal  is  only 
exposed  to  the  risk  of  "observed"  death  for  time  r . If  this  happens 
then  p * r . A usual  assumption  is  that  r = \ , which  we  have  seen 
is  a reasonable  value  of  p . It  is  perhaps  worth  noting  that  at 
the  quadratic  is  positive  so  that  is  always  less  than  the  posterior 

mode,  hence  there  may  be  a systematic  error  in  using  in  under- 

estimating the  death-rate. 


Further  progress  may  be  made  by  returning  to  the  density  of  £ and 
H , Equation  (12),  and  considering  the  second  moments.  Clearly  £ has 
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where  <i>  is  the  solution  of  (14).  (Note:  d + w = 1 .)  The  form  of 
this  variance  is  particularly  interesting.  It  consists  of  two  parts, 
derived  respectively  from  the  variances  of  S and  n , the  first  of 
which  tends  to  zero  as  N 00  , but  the  second  of  which  does  not 
* (unless  the  proportion  of  deaths  or  withdrawals  observed  tends  to  zero). 
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The  latter  part  is  a multiple  of  v , the  prior  variance  of  p . 
Consequently,  however  large  N is  , there  remains  some  uncertainty 
as  to  the  true  death-rate,  ascribable  to  the  original  uncertainty  about 
p . The  estimator  given  by  (14)  is  therefore  not  inconsistent  as  suggested 
by  Breslow  and  Crowley.  Neither  is  it  consistent.  There  is  just  not 
enough  information  in  the  data,  however  much  of  it  there  is,  to  be  sure 
about  the  death-rate.  The  original  uncertainty  about  p is  present 
however  large  is  the  data  set. 

It  is  of  some  interest  to  see  the  effect  of  the  prior  variance  of 
p , v , on  the  very  large  sample  variance  of  <j>  . Letting  N ■>  00  in 
(17),  approximating  $ by  , r by  p and  using  the  facts  that  D/N 
and  S/N  tend  to  the  chances  given  in  (5) , tedious  algebra  establishes 
that  the  variance  of  <p  tends  to 


4>(1  - 4>)8 

(1  - 6)  + p(6  - *) 


2 


v . 


For  example,  if  <p  = 9 so  that  the  withdrawal-  and  death-rates  are 
2 

equal,  this  is  $ v , or  the  standard  deviation  is  scaled  down  from 
U U 

v2  for  p to  (fiv^  for  4>  If,  as  we  suggest  below,  v is  small, 

the  residual  uncertainty  about  <p  is  even  smaller,  especially  for 
low  death-rates. 


r ' 

14 

Discussion: 

f 

Many  statisticians  will  object  to  our  argument  because  it  introduces 
a prior  distribution  for  p , and,  less  importantly,  $ and  0 . 

We  will  discuss  this  in  a moment,  but  first  1 would  ask  them  to  consider 
the  operational  consequences  of  the  method.  Rather  than  committing 
oneself  to  a rigid  assumption  about  F and  H , the  method  allows  for 
flexibility  in  these  functions  and  does  not  involve  a committment  to 
any  particular  form.  The  result  of  this  is  that  instead  of  producing 
estimators  which  are  inconsistent  unless  the  rigid  assumption  is  exactly 
true,  information  is  provided  about  the  uncertainty  in  the  death-rate 
and  the  role  played  in  this  uncertainty  by  the  lack  of  knowledge  of  p . 

Our  method  therefore  seems  both  flexible  and  realistic. 

Here  is  no  place  to  argue  in  favor  of  p having  a distribution. 
Essentially  without  it  inconsistency  of  judgments  cannot  be  avoided. 

The  choice  of  the  exact  form  of  f is  not  a subject  that  a statistician 
should  consider  on  his  own:  he  must  necessarily  consult  with  the 
demographer  or  reliability  engineer  who  is  familiar  with  the  situation. 

The  statistician's  task  is  to  help  the  scientist  to  articulate  his 
knowledge  of  the  problem  in  probabilistic  terms.  Here  it  may  not  be 
unreasonable  to  assume  p has  a fairly  tight  distribution  around  p = \ , 
since,  as  we  saw  above,  many  familiar  situations  give  values  of  p 
either  exactly  h or  very  near  to  it;  and  p near  0 or  1 only 
arises  in  extreme  situations.  But  there  is  one  case  where  p is  near 
3/4  , namely  in  calendar-year  counts  of  academic  staff,  since  most 
withdrawals  take  place  at  the  end  of  the  academic  year  in  September. 

In  our  large-sample  analytic  treatment  only  the  mode  and  variance  of  p 
matter.  With  a uniform  distribution  of  p , v = 1/12  , and  values  sub- 
stantially less  than  this  seem  in  order. 
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We  have  made  the  severe  assumption  that  p is  independent  of 
4>  and  9 . The  appropriateness  of  this  must  be  a matter  for  discussion 
and  cannot  be  decided  on  purely  statistical  grounds.  It  would  not  be 
difficult  to  incorporate  dependence  into  the  argument  but  without  a 
specific  application  in  mind  it  seems  idle  to  speculate  on  the  form 
it  might  take.  In  some  cases  p might  be  a function  of  9 and  $ , 
as  when  F and  H are  both  exponential  (Equation  (7)  above),  identifi- 
ability  is  ensured  and  consistency  obtainable  though  the  likelihood 
is  rather  complicated. 

We  have  considered  only  a single  interval.  For  several  intervals 
1^  , it  is  easy  to  see  the  likelihood  is  the  product  of  factors  like 
(6)  with  data  , W^  , 3^  (S^  ^ + S^)  and  parameters 

^k  * ®k  * Pk  ’ Consequently  the  likelihood  factors.  However  it  seems 
reasonable  that  the  priors  will  not  factor  since  the  death-rates  » 
for  example,  are  likely  to  be  similar  in  adjacent  intervals.  This 
correlation  will  lead  to  smoothing  of  the  rates  as  a function  of  k , 
but  details  remain  to  be  worked  out. 

I am  grateful  to  Professor  Richard  E.  Barlow  who  invited  me  to 
Berkeley,  introduced  me  to  life-table  problems  and  challenged  me  to 
provide  a Bayesian  analysis.  This  paper  is  the  reply  to  the  challenge. 
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